suppressPackageStartupMessages(require(spdep))
suppressPackageStartupMessages(require(sf))
suppressPackageStartupMessages(require(tidycensus))
suppressPackageStartupMessages(require(dplyr))
suppressPackageStartupMessages(require(tmap))
tmap_mode("view")
tract <- tidycensus::get_acs( # Taking data from census
geography = 'tract',
variables = c("Households_" = 'B22010_001', "SNAP_" = 'B22010_002'),
state="MI",
county=c("Oakland", "Macomb", "Wayne"),
geometry = TRUE,
output = 'wide') %>%
st_transform(4326) %>%
mutate(
pct_snap = SNAP_E/Households_E
# pct_snap = ifelse(is.nan(pct_snap), 0, pct_snap)
# pct_snap = ifelse(is.na(pct_snap), 0, pct_snap)
) %>%
filter(
as.vector(sf::st_area(.)) > 0,
st_coordinates(st_centroid(.))[,2] > 33.667,
!is.na(pct_snap),
!is.nan(pct_snap) # Lets remove the NA values
)
suppressMessages(tmap_mode('plot'))
tm_shape(tract) +
tm_polygons(
col='pct_snap',
palette = 'Blues',
border.col = 'white',
lwd=0.3,
style = 'quantile',
n=5,
contrast=c(0,1),
title = "Share of Households (0 to 1)") +
tm_shape(tract %>% st_union()) + # adding a nice border!
tm_polygons(col = NA, alpha = 0, border.col = 'black') +
tm_layout(main.title = "Households on Food Assistance (SNAP)")
tract_nb <- poly2nb(tract,queen=T) %>%
include.self() # This turns it into Gi*! Otherwise, it's just G
tract_W <- nb2listw(tract_nb,
style="W",
zero.policy = TRUE)
tract_local_g <- spdep::localG(
x = tract$pct_snap,
listw = tract_W,
zero.policy = TRUE) %>%
as.vector()
# I made us another function!
make_g_bin <- function(v){
case_when(
v > -1.96 & v < 1.96 ~ "Insignificant",
v >= 1.96 & v < 2.58 ~ "Hotspot (>1.96)",
v >= 2.58 ~ "Hotspot (>2.58)",
v < -1.96 & v > -2.58 ~ 'Coldspot (<-1.96)',
v < -2.58 ~ 'Coldspot (<-2.58)',
TRUE ~ 'Unknown'
) %>% as.factor() %>%
ordered(c(
"Unknown",
"Coldspot (<-2.58)",
"Coldspot (<-1.96)",
"Insignificant",
"Hotspot (>1.96)",
"Hotspot (>2.58)"))
}
g_palette <- c('grey','#0571B0','#92C5DE','cornsilk','#F4A582','#CA0020')
# make_g_bin(v = tract_local_g)
# Apply our colors on the fly to the 'col' column
tract %>% mutate(G=make_g_bin(tract_local_g)) %>%
tm_shape() +
tm_polygons(
col = "G",
palette = g_palette,
border.col = 'grey80',
lwd=0.4,
title = "Local Getis-Ord G") +
tm_layout(main.title = "Households on Food Assistance (SNAP)")
tract %>% mutate(G=make_g_bin(tract_local_g)) %>%
tm_shape() +
tm_polygons(
col = "G",
palette = g_palette,
border.col = 'grey80',
alpha = 0.5,
lwd=0.4,
title = "Local Getis-Ord G") +
tm_shape(st_union(tract)) +
tm_polygons(zindex = 200,lwd=2,alpha = 0,border.col = "black")+
tm_tiles(leaflet::providers$Stamen.TerrainLabels) +
tm_layout(main.title = "Households on Food Assistance (SNAP)")
suppressMessages(tmap_mode("view"))
tmap_options(basemaps = c("Stamen.TonerLite","Stamen.Terrain","Esri.WorldTopoMap","Esri.WorldImagery"))
tract %>% mutate(G=make_g_bin(tract_local_g)) %>%
tm_shape() +
tm_polygons(
col = "G",
palette = g_palette,
border.col = 'grey80',
alpha = 0.5,
lwd=0.4,
title = "Local Getis-Ord G",
popup.vars=c("Tract: " = "NAME", "Percent of households on food assistance (x 100)" = "pct_snap", "Classification: " = "G"))+
tm_shape(st_union(tract)) +
tm_polygons(zindex = 200,lwd=2,alpha = 0,border.col = "black")+
tm_tiles(leaflet::providers$Stamen.TerrainLabels) +
tm_layout(main.title = "Detroit Metro Households on Food Assistance")
Conclusion: Detroit and its surrounding areas were selected for local Getis-Ord G and local Moran’s I spatial analysis. Data of Oakland, Macomb, Wayne counties in the state of Michigan were collected from the census to be analyzed. In our map, we have classified red as hotspots i.e. areas where red is where people use the most food stamps, here the food stamp data is represented as SNAP in the data sheet, and on the other hand, blue is meant as coldspots meaning people in those areas using less stamps for food. We see from the map that there is a large clustering of red in the Detroit area and a small clustering in the Pontiac and Inkster areas. On the other hand, areas like Rochester, Milford, Wixom, Berkley, Clawson, Farmington, Plymouth, Livonia etc. have high and low clustering in blue. That is, most people in the Detroit region depend on food stamps provided by the government to meet their food needs.
tract_nb <- poly2nb(tract,queen=T)
tract_W <- nb2listw(tract_nb,
style="W",
zero.policy = FALSE)
tract_local_i <- spdep::localmoran(
x = tract$pct_snap,
listw = tract_W) %>%
as_tibble()
tract_local_i
## # A tibble: 1,164 × 5
## Ii E.Ii Var.Ii Z.Ii `Pr(z != E(Ii))`
## <localmrn> <localmrn> <localmrn> <localmrn> <localmrn>
## 1 0.3722461 -2.656096e-04 0.051292905 1.6447932 0.1000124552
## 2 -0.1895189 -1.044032e-04 0.017269309 -1.4413706 0.1494800072
## 3 0.2504802 -2.495911e-04 0.036087732 1.3198547 0.1868835331
## 4 1.9432558 -1.302028e-03 0.301675376 3.5403886 0.0003995383
## 5 2.8093032 -7.619053e-03 0.971157312 2.8584465 0.0042572084
## 6 -0.2993126 -5.344158e-04 0.123917467 -0.8487554 0.3960174002
## 7 0.3604603 -3.315454e-04 0.042570536 1.7486472 0.0803520271
## 8 0.3669423 -3.425703e-04 0.079448609 1.3030465 0.1925588791
## 9 1.4714524 -1.059055e-03 0.245438980 2.9722608 0.0029561550
## 10 0.1385828 -6.308377e-06 0.001463526 3.6226720 0.0002915754
## # … with 1,154 more rows
# Transferring values from test into our sf object
tract$lisa_I <- tract_local_i$Ii
tract$lisa_EI <- tract_local_i$E.Ii
tract$lisa_Z <- tract_local_i$Z.Ii
tract$lisa_p <- p.adjust(tract_local_i$`Pr(z != E(Ii))`,
method = "fdr")
# A functino to bin I values
make_i_bin <- function(v){
case_when(
v <= 0.05 ~ "Significant",
v > 0.05 ~ "Insignificant",
is.na(v) ~ "Unknown",
TRUE ~ 'Unknown'
) %>% as.factor() %>%
ordered(c(
"Unknown",
"Insignificant",
"Significant"))
}
require(ggplot2)# load ggplot2 package
## Loading required package: ggplot2
tract$lag_snap <- spdep::lag.listw(x = tract_W, var = tract$pct_snap)
# creating moran_plot function
moran_plot = function(var, w_var, xlab, ylab, main_title = NULL, normalize = T){
if(missing(main_title)){
main_title = "Moran's I Plot"
}
if(normalize){
var = scale(var) %>% as.vector
w_var = scale(w_var) %>% as.vector()
}
# lebel the values as "HH", "HL", "LL", "LH"
plot_dat = data.frame(var,w_var,stringsAsFactors = F) %>%
mutate(quad = case_when(
var > 0 & w_var > 0 ~ "HH",
var > 0 & w_var < 0 ~ "HL",
var < 0 & w_var < 0 ~ "LL",
var < 0 & w_var > 0 ~ "LH",
TRUE ~ 'Unknown'))
#plotting the data in quadrants
plot_dat %>% ggplot() +
geom_vline(xintercept = 0) +
geom_hline(yintercept = 0) +
geom_point(aes(x=var,y=w_var,color=quad)) +
xlab(xlab) +
ylab(ylab) +
geom_smooth(
mapping = aes(x=var,y=w_var),
formula = 'y ~ x',
method = 'lm',
color="black",
size=0.5,
se=F) +
scale_color_manual(
name = "",
values = c('Unknown' = 'grey','HH'="firebrick",'HL'="brown1",'LL'="navy",'LH'="cornflowerblue")) +
ggtitle(main_title)
}
moran_plot( #giving values to the function
var = tract$pct_snap,
w_var = tract$lag_snap,
xlab = "Households on SNAP (%)",
ylab = "Lagged Households on SNAP (%)")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
The Moran’s I plot here suggests that there is a positive spatial autocorrelation between SNAP usage which means that if SNAP usage increases, the SNAP usage of its neighbors will increase too.
tract %>% tm_shape() +
tm_polygons(
col = "lag_snap",
palette = '-RdBu',
border.col = 'grey80',
lwd=0.4,
style = 'cont',
title = "Lagged Households on SNAP (%)") +
tm_layout(main.title = "Households on Food Assistance (SNAP)")
make_i_quads <- function(v, w_v, p, normalize = T){
if(normalize){
v = scale(v) %>% as.vector
w_v = scale(w_v) %>% as.vector()
}
case_when(
p > 0.05 ~ "Insignificant", ### Make sure that 'p' has been adjusted!
v < 0 & w_v < 0 ~ "Low-Low",
v < 0 & w_v > 0 ~ "Low-High",
v > 0 & w_v < 0 ~ "High-Low",
v > 0 & w_v > 0 ~ "High-High",
TRUE ~ 'Unknown'
) %>% as.factor() %>%
ordered(c(
"Insignificant",
"Low-Low",
"Low-High",
"High-Low",
"High-High"))
}
i_palette <- c('cornsilk','#0571B0','#92C5DE','#F4A582','#CA0020')
tract %>% mutate(quad = make_i_quads(v = pct_snap, w_v = lag_snap,p = lisa_p)) %>%
tm_shape() +
tm_polygons(
col = "quad",
palette = i_palette,
border.col = 'grey70',
lwd=0.4,
title = "Local Moran's I") +
tm_layout(main.title = "Local Moran's I: SNAP Benefit Usage in LA County")
Conclusion: From the local Moran-Eye map, we see that there is a greater clustering of food in the Detroit area, and some blue areas to the north and northwest of the red zone. Historically, Detroit was a very successful auto industrial city. But in the mid-1900s, competition in the auto industry increased. In order to compete with Japanese and German auto companies, these companies began to relocate to other cities to build new, more efficient plants. Chrysler, one of Detroit’s dominant auto industries, went bankrupt in 1979. All these economic chaos has eventually made Detro a dying city and currently it has the highest unemployment rate in the entire USA (over 16%). So people here are more dependent on food stamps because they don’t have enough money to buy food. Besides, the wealthier white population migrated north and northwest of Detroit due to the economic crisis, declining property values, lack of jobs, etc., which indicates the reason for the blue clustering of those places.